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Abstract 

We present a numerical investigation of the ray dynamics in a paraxial optical cavity 
when a ray splitting mechanism is present. The cavity is a conventional two- mirror stable 
resonator and the ray splitting is achieved by inserting an optical beam splitter perpendic- 
ular to the cavity axis. We show that depending on the position of the beam splitter the 
optical resonator can become unstable and the ray dynamics displays a positive Lyapunov 
exponent. 

PACS numbers: 42.60.Da, 42.65. Sf, 42.15.-i 

1 INTRODUCTION 

A beam splitter (BS) is an ubiquitous optical device in wave optics experiments, used e.g., 
for optical interference, homodyning, etc. In the context of geometrical optics, light rays are 
split into a transmitted and reflected ray by a BS. Ray splitting provides an useful mechanism 
to generate chaotic dynamics in pseudointegrable^ and soft-chaoticpEI closed systems. In this 
paper we exploit the ray splitting properties of a BS in order to build an open paraxial cavity 
which shows irregular ray dynamics as opposed to the regular dynamics displayed by a paraxial 
cavity when the BS is absent. 

Optical cavities can be classified as stable or unstable depending on the focussing properties 
of the elements that compose it^ An optical cavity formed by 2 concave mirrors of radii R 
separated by a distance L is stable when L < 2R and unstable otherwise. If a light ray is 



injected inside the cavity through one of the mirrors it will remain confined indefinitely inside 
the cavity when the configuration is stable but it will escape after a finite number of bounces 
when the cavity is unstable (this number depends on the degree of instability of the system). 
Both stable and unstable cavities have been extensively investigated since they form the basis 
of laser physics.^ Our interest is in a composite cavity which has both aspects of stability 
and instability. The cavity is made by two identical concave mirrors of radii R separated by 
a distance L, where L < 2R so that the cavity is globally stable. We then introduce a beam 
splitter (BS) inside the cavity, oriented perpendicular to the optical axis (FigQ). In this way 
the BS defines two subcavities. The main idea is that depending on the position of the BS the 
left (right) subcavity becomes unstable for the refiected rays when Li (L2) is bigger than R, 
whereas the cavity as a whole remains always stable {Li + L2 < 2R) (Fig. [2|). 

Our motivation to address this system originates in the nontrivial question whether there 
will be a balance between trapped rays and escaping rays. The trapped rays are those which 
bounce infinitely long in the stable part of the cavity, while the escaping ones are those which 
stay for a finite time, due to the presence of the unstable subcavity. If such balance exists it 
could eventually lead to transient chaos since it is known in literature that instability (positive 
Lyapunov exponents) and mixing (confinement inside the system) form the skeleton of chaos.l^ 
The BS is modelled as a stochastic ray splitting elemenlPby assuming the refiection and trans- 
mission coefficients as random variables. Within the context of wave optics this model cor- 
responds to the neglect of all interference phenomena inside the cavity; this would occur, for 
instance when one injects inside the cavity a wave packet (or cw broad band light) whose lon- 
gitudinal coherence length is very much shorter than the smallest characteristic length of the 
cavity. The stochasticity is implemented by using a Monte Carlo method to determine whether 
the ray is transmitted or refiected by the BS.I^When a ray is incident on the ray splitting surface 
of the BS, it is either transmitted through it with probability p or refiected with probability 
1 — p, where we will assume p = 1/2, i.e. we considered a 50%/50% beam splitter (Fig El)- 
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We then follow a ray and at each reflection we use a random number generator with a uniform 
distribution to randomly decide whether to reflect or transmit the incident ray. 
Our system bears a close connection with the stability of a periodic guide of paraxial lenses as 
studied by LonghiP While in his case a continuous stochastic variable e„ represents a pertur- 
bation of the periodic sequence along which rays are propagated, in our case we have a discrete 
stochastic parameter Pn which represents the response of the BS to an incident ray. As will be 
shown in section II, this stochastic parameter can take only two values, either 1 for transmitted 
rays or -1 for reflected ray; in this sense, our system (as displayed in FigE} is a surprisingly 
simple realization of a bimodal stochastic paraxial lens guide. 

The structure of the paper is as follows. In section II we describe the ray limit, and the 
paraxial map or ABCD matrix associated with rays that propagate very close to the axis of 
the cavity. In section III we present the results of the numerical simulations for the paraxial 
map associated with our ray optical system; these simulations are based on standard numerical 
tools developed in non-linear dynamics theory. Finally, in section IV, we detail the conclusions 
of our work. 

2 Ray Dynamics and the Paraxial Map 

The time evolution of a laser beam inside a cavity can be approximated classically by using 
the ray optics limit, where the wave nature of light is neglected. Generally, in this limit the 
propagation of light in a uniform medium is described by rays which travel in straight lines, 
and which are either sharply reflected or refracted when they hit a medium with a different re- 
fractive index. To fully characterize the trajectory of a ray in a strip resonator or in a resonator 
with rotational symmetry around the optical axis, we choose a reference plane z = constant 
(perpendicular to the optical axis z), so that a ray is specified by two parameters: the height q 
above the optical axis and the angle 6 between the trajectory and the same axis. Therefore we 
can associate a ray of light with a two dimensional vector r = (g, 6). This is illustrated in the 
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two mirror cavity show in Fig. 01 where the reference plane has been chosen to coincide with 
the beam sphtter (BS). Given such a reference plane which is also called Poincare Surface 
of Section (SOS)/ a round trip (evolution between two successive reference planes) of the ray 
inside the cavity can be calculated by the monodromy matrix M„, in other words r^+i = M„r„, 
where the index n determines the number of round trips. The monodromy matrix M„ describes 
the linearized evolution of a ray that deviates from a reference periodic orbit. A periodic orbit 
is said to be stable if |TrM„| < 2. In this case nearby rays oscillate back and forth around the 
stable periodic orbit with bounded displacements both in q and 9. On the other hand when 
|TrM„| > 2 the orbit is said to be unstable and rays that are initially near this reference orbit 
become more and more displaced from it. 

For paraxial trajectories, where the angle of propagation relative to the axis is taken to be 
very small (i.e. sin(0) = tan(0) = 9), the reference periodic trajectory coincides with the optical 
axis and the monodromy matrix is identical to the ABCD matrix of the system. The ABCD 
matrix or paraxial map of an optical system is the simplest model one can use to describe 
the discrete time evolution of a ray in the optical system.* Perhaps the most interesting and 
important application of ray matrices comes in the analysis of periodic focusing (PF) systems 
in which the same sequence of elements is periodically repeated many times down in cascade. 
An optical cavity provides a simple way of recreating a PF system, since we can think of a 
cavity as a periodic series of lenses (see Fig|3]). In the framework of geometric ray optics, PF 
systems are classified, optical cavities, as either stable or unstable. 

Without essential loss of generality we restrict ourselves to the case of a symmetric cavity (i.e. 
two identical spherical mirrors of radius of curvature R). We take the SOS coincident with the 
surface of the BS. After intersecting a given reference plane z,, a ray is transmitted (reflected), 
it will undergo a free propagation over a distance L2 (Li), followed by a reflection on the curved 
mirror M2 (Mi), and continue propagating over the distance L2 (Li), to hit the surface of the 
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beam splitter again at Zj+i. In Fig|3]the sequence of Zi represents the successive reference planes 
after a round trip. In the paraxial approximation each round trip (time evolution between two 
successive intersections of a ray with the beam splitter) is represented by: 



6'n+l = Cqn + Drfir, 



(1) 



where 



A„ = 1 - 2LrjR, 




I-LJR), 



C 



2/R, 



Dn = l- 2L JR 



and 



L 



n 



L + Pna 
2 



We have defined L = Li + L2 and a = L2 — Li, the stochastic parameter pn = ±1 determines 
whether the ray is transmitted (p„ = 1) or is reflected (p„ = —1). 

The elements of the ABCD matrix depend on n because of the stochastic response of the 
BS, which determines the propagation for the ray in subcavities of different length (either Li 
or L2). In this way a random sequence of reflections {pn = 1) and transmissions (p„ = — 1) 
represents a particular geometrical realization of a focusing system. If we want to study the 
evolution of a set of rays injected in the cavity with different initial conditions {qo,0o), we 
have two possibilities, either use the same random sequence of reflections and transmissions 
for all rays in the set or use a different random sequence for each ray. In the latter case, we 
are basically doing an ensemble average over different geometrical configurations of focusing 
systems. As we shall see later it is convenient, for computational reasons, to adopt the second 
method. 

In the next section we report several dynamical quantities that we have numerically calcu- 
lated for paraxial rays in this system, using the map described above (Eq^ . The behavior 



5 



of these quantities, namely, the SOSs, the exit basins, the Lyapunov exponent and the escape 
rate, is analyzed as a function of the displacement (A) of the BS with respect to the center of 
the cavity (see FigH]). 

3 Results 

The paraxial map of Eq^ describes an unbounded system, that is rays are allowed to go in- 
finitely far from the cavity axis. In order to describe a physical paraxial cavity we have to 
keep the phase space bounded, i.e. it is necessary to artificially introduce boundaries for the 
position and the angle of the ray.l^ The phase space boundaries that we have adopted to decide 
whether a ray has escaped after a number of bounces or not is the beam waist (wq) and the 
diffraction half-angle (Bq) of a gaussian beam confined in a globally stable two-mirror cavity. 
Measured at the center of the cavity, Wq = —^^^^\f^^^ and the corresponding diffraction 
half-angle Bq = arctan(^^^).'^ For our cavity configuration we assume R = 0.15m , L = 0.2m 
and Xiight = 500nm, from which follows that wq = 5.3 x 10~^m and Bq = 0.15 x lO^^rad. 
One should keep in mind that this choice is somewhat arbitrary and other choices are certainly 
possible. The effect of this arbitrariness on our results will be discussed in detail in section D. 

3.1 Poincare surface of section (SOS) 

We have first calculated the SOS for different positions of the BS. In order to get a qualitative 
idea of the type of motion, we have chosen as transverse phase space variables y = q and 
Vy = sin{6) ^ 9. The successive intersections of a trajectory with initial transverse coordinates 
go = 1 X lO^^m, = are represented by the different black points in the surface of section. 
The different SOSs are shown in Fig El In Fig0 (a) we show the SOS for A = 0, while in (b) 
A = 1 X 10^'^m and in (c) A = 2 x lO^^m. In (a) it is clear that the motion is completely 
regular (non-hyperbolic); the on-axis trajectory represents an elliptic fixed point for the map. 
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In (b), where the BS is shghtly displaced from the center (A = 1 x lO^^m) we can see that this 
same trajectory becomes unstable because of the presence of the BS, and spreads over a finite 
region of the phase space to escape after a large number of bounces (n = 5 x 10^). In this case 
we may qualify the motion as azimuthally ergodic. The fact that the ray-splitting mechanism 
introduced by the BS produces ergodicity is a well known resultP for a closed billard. We find 
here an analogue phenomenon, with the difference that in our case the trajectory does not ex- 
plore uniformly (but only azimuthally) the available phase space, because the system is open. 
Finally, in (c) we see that the fixed point in the origin becomes hyperbolic, and the initial orbit 
escapes after relatively few bounces (n = 165). 

3.2 Exit basin diagrams 

It is well known that chaotic hamiltonian systems with more than one exit channel exhibit 
irregular escape dynamics which can be displayed, e.g., by plotting the exit basin.l^ For our 
open system we have calculated the exit basin diagrams for three different positions of the 
BS (Figini). These diagrams can be constructed by defining a fine grid (2200 x 2200) of initial 
conditions (go, ^o)- We then follow each ray for a sufficient number of bounces so that it escapes 
from the cavity. When it escapes from above (0„ > 0) we plot a black dot in the corresponding 
initial condition, whereas when it escapes from below {6n < 0) we plot a white dot. 
In FiglBl (a) we show the exit basins for A = 0.025m, the uniformly black or white regions of 
the plot correspond to rays which display a regular dynamics before escaping, and the dusty 
region represents the portion of phase space where there is sensitivity to initial conditions. In 
Fig. ini(b), we show the same plot for A = 0.05m, and in (c) for A = 0.075m. 
The exit basins plots in FigEl illustrate how the scattering becomes more irregular as the BS is 
displaced from the center. In particular, we see how regions of regular and irregular dynamics 
become more and more interwoven as A increases. Instead, for small values of A as in FiglHfa), 
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we can see that there is a single dusty region with a uniform distribution of white and black 
dots in which no islands of regularity are present. 

3.3 Escape rate and Lyapunov exponent 

The next dynamical quantities we have calculated are the escape rate 7 and the Lyapunov 
exponent A. The escape rate is a quantity that can be used to measure the degree of openness 
of a system.^ For hard chaotic systems (hyperbolic), the number of orbits still contained in 
the phase space after a long time (measured in number of bounces n) decreases as A^o exp(— 7n), 
while for soft chaotic systems, the stickiness to Kolmogorov-Arnold-Moser (KAM) islands (or 
islands of stability) leads to a power law decay Agn"'^.'^ The Lyapunov exponent is the rate 
of exponential divergence of nearby trajectories. 

Since both A and 7 are asymptotic quantities they should be calculated for very long times. 
In our system long living trajectories are rare, and in order to pick them among the grid of 
initial conditions A^o one has to increase Aq beyond the computational capability. To overcome 
this difficulty we choose a different random sequence for each initial condition. In this way 
we greatly increase the probability of picking long living orbits given by particularly stable 
random sequences. These long living orbits in turn make possible the calculation of asymptotic 
quantities such as A or 7. 

The escape rate 7 was determined measuring A^„, as the slope of a linear fit in the Nn/No 
versus n curve, in a logarithmic scale; the total number of initial conditions A"o being chosen 
as 2200 X 2200. 

We have calculated the dependence of 7 with the displacement A of the BS from the center of 
the cavity, where < A < L/2. Since for A > R — L/2 the left subcavity becomes unstable, 
it would seem natural to expect that this position of the BS would correspond to a critical 
point. However, we have found by exphcit calculation of both the Lyapunov exponent and 
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the escape rate, that such a critical point does not manifest itself in a sharp way, rather we 
have observed a finite transition region (as opposed to a single point) in which the functional 
dependence of A and 7 change in a smooth way. In Fig U\ (a) we show the typical behavior 
of ^vs n in semi-logarithmic plot for three different positions of the BS. The displacement 
of the BS is A = 0.0875m, 0.05m and 0.03125m, and the corresponding slopes (escape rate 7 
measured in units of the inverse number of bounces n) of the linear fit are 7 = 0.17693?2~^, 
0.05371n~^ and 0.01206/2^^ respectively. We have found that the decay is exponential only up 
to a certain time (approximately 70-1000 bounces depending on the geometry of the cavity) 
due the discrete nature of the grid of initial conditions. 

In Fig[7| (b) we see that 7 increases with A, revealing that for more unstable configurations 
there is a higher escape rate, as expected. Its also interesting to notice that the exponential 
decay fits better when the beam splitter is further from the center position, since this leads to 
smaller stability of the periodic orbits of the system. However, the dependence of the escape 
rate with the position of the BS is smooth and reveals that the only critical displacement, where 
the escape rate becomes positive, is A = 0. 

As a next step, we have calculated the Lyapunov exponent A for the paraxial map; A 
is a quantity that measures the degree of stability of the reference periodic orbit. For a two- 
dimensional hamiltonian map there are two Lyapunov exponents (Ai, A2) such that A1 + A2 = 0. 
In the rest of the paper we shall indicate with A the positive Lyapunov exponent which quantifies 
the exponential sensitivity to the initial conditions. We have calculated A for the periodic orbit 
on axis, using the standard techniques,'^ and we have found that the Lyapunov exponent grows 
from zero with the distance of the BS to the center (Fig[3(c)). Therefore, the only critical point 
revealed by the ray dynamics is again the center of the cavity (A = 0), where the magnitudes 
change from zero to a positive value. This result also shows that the presence of the BS with its 
stochastic nature introduces exponential sensitivity to initial conditions in the system for every 
A 7^ 0, even when both subcavities are stable. This surprising fact can be explained by taking 
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into account the well known probabilistic theorem by Furstenberg on the asymptotic limit of 
the rate of growth of a product of random matrices (PRM).'^From this theorem we expect that 
the asymptotic behavior of the product M„ of a uniform sequence uj of independent, random, 
unimodular, D x D matrices, and for any nonzero vector y G 9?^: 

lim -(ln|M„y|) = Ai > 0, (2) 



where Ai is the maximum Lyapunov characteristic exponent of the system, and the angular 
bracket indicates the average over the ensemble Q of all possible sequences uj. This means that 
for PRM the Lyapunov exponent is a nonrandom positive quantity. In general, it can be said 
that there is a subspace Q* of random sequences which has a full measure (probability 1) over 
the whole space of sequences Q for which nearby trajectories deviate exponentially at a rate Ai. 
Although there exist very improbable sequences in Q which lead to a different asymptotic limit, 
they do not change the logarithmic average (Eql2])P^ We have verified this result, calculating 
the value of A for different random sequences Ui, in the asymptotic limit n = 100000 bounces, 
and we obtained in all cases the same Lyapunov exponent. 

3.4 Mixing properties 

Dynamical randomness is characterized by a positive Kolmogorov- Sinai (KS) entropy per unit 
time hxs^ In closed systems, it is known that dynamical randomness is a direct consequence 
of the exponential sensitivity to initial conditions given by a positive Lyapunov exponent. On 
the other hand, in open dynamical systems with a single Lyapunov exponent A, the exponential 
sensitivity to initial conditions can be related to hxs through the escape rate 7, by the relation!^ 

A = hxs + 7- (3) 

This formula reveals the fact that in an open dynamical system the exponential sensitivity to 
initial conditions induces two effects: one is the escape of trajectories out of the neighborhood of 
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the unstable reference periodic orbit at an exponential rate 7, and the other one is a dynamical 
randomness because of transient chaotic motion near this unstable orbit This dynamical 
randomness is a measure of the degree of mixing of the system and as mentioned before is 
quantified by h^s- Therefore, for a given A, the larger the mixing is, the smaller the escape 
rate, and vice versa. From FigsIZ|^b,c) it is evident that the Lyapunov exponent and the escape 
rate have the same smooth dependence on the BS displacement A and that 7 < A. We have 
calculated the difference A — 7 > for our system and the result is shown in Fi^ (d). 
The actual value of 7(A) depends, for a fixed value of A, on the size of the phase space 
accessible to the system,!^ that is, it depends on wq and Oq. We verified this behavior by 
successively decreasing wq and 6q by factors of 10 (see Table H}, and calculating 7 for each 
of these phase space boundaries. It is clear from these results that 7 increases when the size 
of phase space decreases; in fact for wq, 9q ^ 0, one should get A ~ 7 and the cavity mixing 
property should disappear. It is important to notice that the increment of 7 with the inverse 
of the size of the accessible phase space is a general tendency, independent from the arbitrarily 
chosen boundaries. 





xlO° 




xlO-2 


xlO-3 
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0.17639 


0.17596 


0.19559 


0.25259 



Table 1: Escape rate for different phase space boundaries. As the boundary shrinks 7(A) 
tends to the corresponding value of A(A) = 0.29178n~^. In these calculations the displacement 
of the BS was A = 0.0875m. 

It is important to stress that, although the randomness introduced by the stochastic BS is 
obviously independent from the cavity characteristics, A and 7 show a clear dependance on the 
BS position. When the BS is located at the center of the cavity it is evident for geometrical 
reasons that the ray splitting mechanism becomes ineffective: A = = 7. These results confirm 
what we have already shown in the SOS (Fig. 4). 
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4 CONCLUSIONS 

We have been able to characterize the ray dynamics of our optical cavity with ray splitting by 
using standard techniques in nonlinear dynamics. In particular we have found, both through 
the SOS and the exit basin diagrams, that the stochastic ray splitting mechanism destroys the 
regular motion of rays in the globally stable cavity. The irregular dynamics introduced by the 
beam splitter was quantified by calculating the Lyapunov exponent A; it grows from zero as 
the beam splitter is displaced from the center of the cavity. Therefore, the center of the cavity 
constitutes the only point where the dynamics of the rays is not affected by the stochasticity of 
the BS. The escape rate 7 has been calculated and it has revealed a similar dependence with 
the position of the beam splitter to that of A. Furthermore, we have verified that the absolute 
value of the escape rate tends to that of the Lyapunov exponent as the size of the available 
phase space goes to zero. This result confirms the fact that the escape rate and therefore the 
mixing properties of a map depend sensitively on the choice of the boundary.!^ Because of this 
dependence we cannot claim that our system is chaotic, despite the positiveness of A. However, 
in a future publication we shall demonstrate that ray chaos can be achieved for the same class 
of optical cavities when non-paraxial ray dynamics is allowed.^ 

This project is part of the program of FOM and is also supported by the EU under the IST- 
ATESIT contract. We thank S. Oemrawsingh for useful contributions to software development. 
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Figure 1: Schematic diagram of the cavity model. Two subcavities of length Li and L2 are 
coupled by a BS. The total cavity is globally stable for L = Li + L2 < 2R. A = Li — L/2 
represents the displacement of the BS with respect to the center of the cavity. 




(c) 




Figure 2: The different positions of the beam splitter determine the nature of the subcavities. 
In (a) the BS is in the middle, so the 2 subcavities are stable, in (b) the left cavity is unstable 
and the right one is stable, and (c) the unstable (stable) cavity is on the left (right) (b). 
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Figure 3: A ray on a reference plane [z = const) perpendicular to the optical axis (Z) is 
specified by two parameters: the height q above the optical axis and the angle 6 between the 
direction of propagation and the same axis. When a ray hits the surface of the BS, which we 
choose to coincide with the reference plane, it can be either reflected or transmitted with equal 
probability. For a 50%/50% beam splitter p = 1/2. 



Z2 
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Figure 4: A ray bouncing inside an optical cavity can be represented by a sequence of lenses 
of focus f = 2/R, followed by a free propagation over a distances L„. Due to the presence of 
the BS, the distance L„ varies stochastically between Li or L2. 
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Figure 5: SOS for (a) A = the ray does not escape, (b) A = 0.001, the ray escapes after 
n = 5 X 10^ bounces and (c) A = 0.02, the ray escapes after n = 165. 
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Figure 7: (a) Linear fits used to calculate the escape rate for three different geometrical 
configurations of the cavity given by A = 0.03125m, A = 0.05m and A = 0.0875m. The time is 
measured in number of bounces [n]. The slope 7 is in units of the inverse of time [n~^]. Fig (b) 
shows the escape rate 7 [n~^] as a function of A. Fig. 
exponents A [n~^] as the BS moves from the center A 
A = 0.10m. Fig. (d) shows the difference between A - 
function. 
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